!>\brief
!!
!! Linear inertia-gravity wave test case.
!!
!! \n
!!
!! This is a 2D test case described in <a
!! href="http://dx.doi.org/10.1175/1520-0493(1994)122<2623:EAAOTK>2.0.CO;2">[Skamarock,
!! Klemp, Mon. Wea. Rev., 1994]</a>, <a
!! href="http://dx.doi.org/10.1002/fld.1392">[Ahmad, Lindeman, Int. J.
!! Numer.  Methods Fl., 2007]</a>, <a
!! href="http://dx.doi.org/10.1016/j.jcp.2007.12.009">[Giraldo,
!! Restelli, J. Comp. Phys., 2008]</a> and <a
!! href="http://dx.doi.org/10.1137/070708470">[Restelli, Giraldo, SIAM
!! J. Sci. Comp., 2009]</a>.
!!
!! \note This test uses periodic boundary conditions in the \f$x\f$
!! direction.
!<----------------------------------------------------------------------
module mod_gravity_wave_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_gravity_wave_test_constructor, &
   mod_gravity_wave_test_destructor,  &
   mod_gravity_wave_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof, t_ref,  &
   coeff_neu,        &
   coeff_visc,       &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

 ! public interface
 character(len=*), parameter   :: test_name = "gravity-wave"
 character(len=100), protected :: test_description(5)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 character(len=*), parameter   :: ref_prof = "isothermal"
 real(wp), target              :: t_ref

 logical, protected :: mod_gravity_wave_test_initialized = .false.
 ! private members
 real(wp) :: &
   nn,     & !< uniform Brunt-V&auml;is&auml;l&auml; frequency
   t_s,    & !< surface temperature
   dtheta, & !< temperature perturbation
   xc,     & !< perturbation position
   hh,     & !< domain vertical size
   aa,     & !< perturbation width
   u_mean    !< uniform flow
 character(len=*), parameter :: &
   test_input_file_name = 'gravity_wave.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_gravity_wave_test'

 ! Input namelist
 namelist /input/ &
   nn, t_s, dtheta, xc, hh, aa, u_mean 

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_gravity_wave_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_gravity_wave_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Inertia-gravity wave test, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  disturbance dtheta = ',dtheta,';'
   write(test_description(3),'(a,e9.3,a)') &
     '  mean flow velocity = ',u_mean,';'
   write(test_description(4),'(a,e9.3,a)') &
     '  Brunt-Vaisala frequency = ',nn,';'
   write(test_description(5),'(a,e9.3,a)') &
     '  surface temperature = ',t_s,'.'

   ! All the physical constants use the default values: phc unchanged

   t_ref = 300.0_wp

   mod_gravity_wave_test_initialized = .true.
 end subroutine mod_gravity_wave_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_gravity_wave_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_gravity_wave_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_gravity_wave_test_initialized = .false.
 end subroutine mod_gravity_wave_test_destructor

!-----------------------------------------------------------------------

  !> Viscosity
  pure function coeff_visc(x,rho,p,t) result(mu)
   real(wp), intent(in) :: x(:,:), rho(:), p(:), t(:)
   real(wp) :: mu(2,size(x,2))

    mu = 0.0_wp
  end function coeff_visc

!-----------------------------------------------------------------------

 !> Potential temperature perturbation
 !!
 !! The initial disturbance is specified as
 !! \f{displaymath}{ \displaystyle
 !!  \Delta\theta = \Delta\theta_0 \frac{\sin(\pi \frac{z}{H})}
 !!  {1+\left(\frac{x-x_c}{a}\right)^2}.
 !! \f}
 !! where the parameters \f$\Delta\theta_0\f$, \f$H\f$, \f$x_c\f$ and
 !! \f$a\f$ are specified by the input file. Notice that no pressure
 !! correction is applied to the undisturbed hydrostatic profile.
 !<
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  real(wp), parameter :: pi_trig = 3.141592653589793_wp
  real(wp) :: delta
  real(wp), dimension(size(x,2)) :: z, zexp, theta, pi, p, t

   ! first define a stable, constant stability atmosphere
   z = x(size(x,1),:) ! vertical coordinate
   delta = 1.0_wp - phc%gravity**2/(phc%cp*nn**2*t_s)
   zexp = exp( nn**2/phc%gravity * z )
   theta = t_s * zexp
   pi = (1.0_wp-delta)/zexp + delta

   ! potential temperature perturbation
   theta = theta + dtheta * sin(pi_trig*z/hh) &
           / (1.0_wp + ( (x(1,:)-xc)/aa )**2)
   ! no pressure perturbations

   ! translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                                        ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + 0.5_wp*u_mean**2 + phc%gravity*z) ! e
   u0(3,:) = u0(1,:)*u_mean                                        ! Ux
   u0(4:,:) = 0.0_wp                                               ! Uyz
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_init
 
!-----------------------------------------------------------------------

 !> Similar to \c coeff_init, without perturbation
 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  real(wp) :: delta
  real(wp), dimension(size(x,2)) :: z, zexp, theta, pi, p, t

   ! first define a stable, constant stability atmosphere
   z = x(size(x,1),:) ! vertical coordinate
   delta = 1.0_wp - phc%gravity**2/(phc%cp*nn**2*t_s)
   zexp = exp( nn**2/phc%gravity * z )
   theta = t_s * zexp
   pi = (1.0_wp-delta)/zexp + delta

   ! translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                                        ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + 0.5_wp*u_mean**2 + phc%gravity*z) ! e
   u0(3,:) = u0(1,:)*u_mean                                        ! Ux
   u0(4:,:) = 0.0_wp                                               ! Uyz
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_outref

!-----------------------------------------------------------------------

 pure function coeff_neu(x,u_r,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: h(size(x,1),2+size(x,1),size(x,2))

   h = 0.0_wp

 end function coeff_neu

!-----------------------------------------------------------------------

end module mod_gravity_wave_test

